\documentclass[reqno]{amsart}
%\usepackage{txfonts}
\usepackage{CJK}
\usepackage{cite}
\usepackage{mathrsfs}
\usepackage{bm}
\usepackage{graphicx}
\usepackage{epsfig}
\usepackage{epstopdf}
\usepackage{amsmath, amssymb, amsthm, amsmath, mathtools}
\usepackage{enumerate}
\usepackage{hyperref}

\textheight 20.0cm \textwidth 14.0cm
\numberwithin{equation}{section}

\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{claim}[theorem]{Claim}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{example}[theorem]{Example}
\numberwithin{equation}{section}
\allowdisplaybreaks

\makeatletter
\@namedef{subjclassname@2020}{%
	\textup{2020} Mathematics Subject Classification}
\makeatother

\renewcommand{\arraystretch}{1.5}

\newcommand{\differential}{\mathop{}\!\mathrm{d}}
\DeclareMathOperator{\Tail}{Tail}
\newcommand{\relphantom}[1]{\mathrel{\phantom{#1}}\negmedspace {}}
\newcommand{\phantomeq}{\relphantom{=}}
\newcommand{\realset}{\mathbb{R}}
\newcommand{\naturalset}{\mathbb{N}}
\DeclareMathOperator*{\esssup}{ess\,sup}
\DeclareMathOperator*{\essinf}{ess\,inf}
\DeclareMathOperator{\supp}{supp}
\DeclareMathOperator{\divergence}{div}

\def\Xint#1{\mathchoice
    {\XXint\displaystyle\textstyle{#1}}%
    {\XXint\textstyle\scriptstyle{#1}}%
    {\XXint\scriptstyle\scriptscriptstyle{#1}}%
    {\XXint\scriptscriptstyle\scriptscriptstyle{#1}}%
\!\int}
\def\XXint#1#2#3{{\setbox0=\hbox{$#1{#2#3}{\int}$}
    \vcenter{\hbox{$#2#3$}}\kern-.5\wd0}}
\def\ddashint{\Xint=}
\def\dashint{\Xint-}
\newcommand{\bbint}[1]{\dashint_{#1}}


\usepackage[sort&compress, capitalise]{cleveref}

\begin{document}
\title{Nonhomogeneous Nonlocal Weak Harnack Inequality}
\keywords{Nonlocal equations; Weak Harnack estimates; General growth}
\begin{abstract}
     In this paper, we prove two weak Harnack inequalities for weak solutions to nonhomogeneous nonlocal equations under general growth condition. Our approach mainly relies on the spread of pointwise positivity in the spirit of De Giorgi classes. Additionally, a refined energy estimate is employed several times as a key element.
\end{abstract}
\maketitle

\section{Introduction}\label{section: introduction}
Let $\Omega\subset\realset^{N}(N\geq2)$ be a bounded domain. In this paper, we consider the following class of nonlocal equations with general growth
\begin{equation}\label{luf}
    \mathcal{L}u=f \quad\text{in}\ \Omega.
\end{equation}
Here, $\mathcal{L}$ is a nonlocal operator of the form
\begin{equation*}
    \mathcal{L}u(x)\coloneq\mathrm{P.V.}\int_{\realset^{N}}g\left(\frac{|u(x)-u(y)|}{|x-y|^{s}}\right)\frac{u(x)-u(y)}{|u(x)-u(y)|}\frac{K(x,y)}{|x-y|^{s}}\differential y,
\end{equation*}
where $\mathrm{P.V.}$ stands for ``in the principal value sense" and $s\in(0,1)$. $K\colon\realset^{N}\times\realset^{N}\rightarrow (0,\infty]$ is a measurable kernel function fulfilling 
\begin{equation*}
    \frac{\lambda}{|x-y|^{N}}\leq K(x,y)=K(y,x)\leq \frac{\Lambda}{|x-y|^{N}}
\end{equation*}
for some $\Lambda\geq\lambda>0$. 
Moreover, the continuous function $g\colon[0,\infty)\rightarrow[0,\infty)$ is strictly increasing satisfying $g(0)=0$, $\lim_{t\rightarrow\infty}g(t)=\infty$ and 
\begin{equation}\label{G}
    1< p\leq \frac{tg(t)}{G(t)}\leq q<\infty \quad\text{with}\ G(t)=\int_{0}^{t}{g(\tau)}\mathrm{d}\tau,
\end{equation}
where the conditions of $G$ will be introduced in Section \ref{section: preliminaries} in detail.
The nonhomogeneous term $f$ satisfies 
\begin{equation*}
|f(x,u)|\leq d_{1}+d_{2}g(|u|)  \quad\text{for a.a. }x\in\Omega
\end{equation*}
with fixed constants $d_{1}$ and $d_{2}$.

Equation \eqref{luf} is regarded as the nonlocal nonhomogeneous analogue corresponding to the $G$-Laplace equation, which is typically defined as
\begin{equation*}
    -\divergence\Bigl(g(|\nabla u|)\frac{\nabla u}{|\nabla u|}\Bigr)=0.
\end{equation*}
The $G$-Laplace equation has been widely studied in past years. H\"older regularity of solutions to this equation involving measures is derived by Zheng, Feng and Zhang in \cite{Glaplacemeasures2015}. We refer readers to \cite{P. Baroni2015,L. Diening2009,Regularityof2017}. The regularity theory for nonlocal elliptic and parabolic equations has been witnessed remarkable advancements. In the special case $g(t)=t^{p-1}$ and $K=|x-y|^{-N}$, the operator $\mathcal{L}$ in \eqref{luf} can also be represented by $(-\Delta)^{s}_{p}$, which is called fractional $p$-Laplacian, see \cite{byun2024holder,kassmann2011new} for instance. Regarding the regularity theory applicable to this class of problems, Cozzi \cite{2016Regularity} develops local boundedness, Hölder continuous, and a Harnack inequality for weak solutions, utilizing a nonlocal nonhomogeneous energy estimate. The local behavior of weak solutions, including boundedness and H\"older continuity, is investigated by Di Castro et al. \cite{di2016local} in the vein of the De Giorgi-Nash-Moser iteration. Liao in \cite{liao2024nonlocalweakharnackestimates} derives weak Harnack estimates for the nonlocal $p$-Laplace equation. 
For nonlocal parabolic equations with bounded measurable coefficients, Kassmann and Weidner in \cite{kassmann2024parabolicharnackinequalitynonlocal} establish a Harnack inequality for weak solutions, where tail terms appear on both sides of the inequality. Moreover, they focus on solutions to integro-differential equations governed by nonlocal operators characterized by nonsymmetric forms, establishing H\"older regularity and Harnack inequalities in \cite{kassmann2022nonlocal} and \cite{kassmann2024nonlocal}, respectively.

In terms of the higher regularity theory, Brasco et al. in \cite{WOS:000447961200017} provide an explicit H\"older exponent for solutions of the $p$-Laplacian fractional equation in the superquadratic case, while the subquadratic case is studied by Garain-Lindgren in \cite{garain2024higher}. In \cite{iannizzotto2016global}, $C^{\alpha}$-regularity up to boundary to weak solutions of fractional $p$-Laplacian equation is investigated using barrier arguments rather than De Giorgi-Nash-Moser techniques by Iannizzotto, Mosconi and Squassina. Very recently, \cite{bögelein2024gradientregularityspharmonicfunctions,bogelein2024regularity} show that the $(s,p)$-harmonic functions are H\"older continuous to an arbitrary H\"older exponent in (0,1) with the restriction of order $s$ and have fractional differentiability of the gradient. 

When $g(\cdot)$ has a general structure, Fang and Zhang \cite{2022Harnack} employ the expansion of positivity and energy estimates to establish a Harnack inequality; Byun et al. \cite{byun2023nonlocal} give a more simplified proof to obtain Harnack inequality, under no assumptions on $G$ other than \eqref{G}. Besides, \cite{byun2021local} proves the local H\"older regularity without boundary data and any priori bounded assumption; see also \cite{bonder2023global,kim2023supersolutionssuperharmonicfunctionsnonlocal}. As for double phase, \cite{BYUN2022110} shows that weak solutions of nonlocal double phase problems are locally bounded and Hölder continuous. We also refer readers to \cite{buryachenko2020local,de2019holder} and references therein.

In this paper, we aim to establish weak Harnack inequalities to nonhomogeneous nonlocal equations with general growth. Inspired by \cite{2022Harnack}, we use the method known as the expansion of positivity, which tells the propagation of pointwise positivity based on measure information. It's necessary to mention that Kim and Lee \cite{kim2023supersolutionssuperharmonicfunctionsnonlocal} have proved a result similar to Remark \ref{theorem: 1}. Their approach is based on the Moser iteration, while ours is mainly based on an energy estimate that considers the influence of nonhomogeneous term. Different from the typical weak Harnack inequality, which displays the nonlocal tail only on the one side, our result takes into account the feature of nonlocal terms on the positive and negative parts of the weak solutions at both sides of the inequality, respectively.

We introduce the tail-space as follows
\begin{equation*}
    L^{g}_{s}(\realset^{N})=\Bigl\{u\text{ is a measurable function in }\realset^{N}\colon \int_{\realset^{N}}g\Bigl(\frac{|u(x)|}{|x-x_{0}|^{s}}\Bigr)\frac{\differential x}{|x-x_{0}|^{N+s}}<\infty\Bigr\}.
\end{equation*}
Besides, the nonlocal tail of $u$ is defined as
$$\mathrm{Tail}(u;B_{R}(x_{0}))=\int_{\realset^{N}\setminus B_{R}(x_{0})}g\Bigl(\frac{|u(x)|}{|x-x_{0}|^{s}}\Bigr)\frac{\differential x}{|x-x_{0}|^{N+s}}.$$
It is obtained from \cite{byun2021local} that $u\in L^{g}_{s}(\realset^{N})$ if and only if $\mathrm{Tail}(u;B_{R}(x_{0}))$ is finite for any $x_{0}\in\realset^{N}$ and $R>0$.

Now, we proceed to state our main results. As shown in the proof of Lemma \ref{lemma3}, the density of the shrinking set can be estimated by the tail of $u_{+}$. Hence, one can allow a nonlocal term for right-hand side. Besides, we assume $u$ nonnegative in $B_{R}(x_{0})$, and thus the tail of $u_{-}$ needs to be added on the left-hand side.
\begin{theorem}\label{theorem: 2}
    Let $k\geq 0$, $0\leq R\leq 1$ and suppose that $u\in\mathbb{W}^{s,G}(\Omega)$, satisfying $u\geq 0$ in $B_{R}(x_{0})\subset\Omega$, is a weak supersolution to \cref{luf}. Then there exists a constant $\eta\in(0,1)$ depending only on ${s,N,p,q,\lambda,\Lambda,d_{2}}$, such that
    \begin{equation*}
        \begin{aligned}
            &\phantomeq \essinf_{B_{\frac{r}{2}}(x_{0})}u+r^{s}g^{-1}\left(r^{s}\Tail\left(u_{-};B_{R}\left(x_{0}\right)\right)\right)+r^{s}G^{-1}G^{*}\left(d_{1}r^{s}\right)
            \\ & \geq \eta r^{s}g^{-1}\left(r^{s}\Tail\left(u_{+};B_{r}\left(x_{0}\right)\right)\right),
        \end{aligned}
    \end{equation*}
    where $B_{2r}(x_{0})\subset B_{R}(x_{0})$. 
\end{theorem}

The second one is a weak Harnack estimate that takes the nonhomogeneous term into account, stated as followed.
\begin{remark}\label{theorem: 1}
    Let $k\geq 0$, $0\leq R\leq 1$ and suppose that $u\in\mathbb{W}^{s,G}(\Omega)$, satisfying $u\geq 0$ in $B_{R}(x_{0})\subset\Omega$, is a weak supersolution to \cref{luf}. Then there exists a constant $\eta\in(0,1)$ depending only on ${s,N,p,q,\lambda,\Lambda,d_{2}}$, such that
    \begin{equation*}
        \begin{aligned}
            &\phantomeq \essinf_{B_{\frac{r}{2}}(x_{0})}u+r^{s}g^{-1}\left(r^{s}\Tail(u_{-};B_{R}(x_{0}))\right)+r^{s}G^{-1}G^{*}(d_{1}r^{s})
            \\ & \geq \eta r^{s}g^{-1}\left(\bbint{B_{r}(x_{0})}g\left(\frac{u(x)}{r^{s}}\right)\differential x\right),
        \end{aligned}
    \end{equation*}
    where $B_{2r}(x_{0})\subset B_{R}(x_{0})$.
\end{remark}

The article is organized as follows. In Section \ref{section: preliminaries}, several basic concepts, inequalities, and function spaces are introduced. In Section \ref{section: energy estimate}, we establish an energy estimate that makes an effort to obtain the expansion of positivity and density lemmas in Section \ref{section: density lemma}. Finally, we prove the weak Harnack inequalities.

\section{Preliminaries}\label{section: preliminaries}
In this section, we shall introduce the definition of weak solutions, function spaces related to solutions, and give some basic inequalities to be used later.

Let $B_{r}(x_{0})\coloneq\{x\in\realset^{N}\colon|x-x_{0}|<r\}$ be an open ball with center $x_{0}$ and radius $r>0$, and center $x_{0}$ will be omitted when there is no ambiguity. A measurable function  $G\colon[0,\infty)\rightarrow [0,\infty)$ is said to be an $N$-function if it is convex and increasing, and satisfies 
\begin{equation*}
    G(0)=0,\quad \lim_{t\to 0^{+}}\frac{G(t)}{t}=0\quad \text{and}\quad \lim_{t\to\infty}\frac{G(t)}{t}=\infty.
\end{equation*}
For an $N$-function $G\colon[0,\infty)\to [0,\infty)$, 
whose conjugate function $G^{*}\colon[0,\infty)\to [0,\infty)$ is defined by
\begin{equation*}
    G^{*}(t)=\sup_{\tau\geq0}\{\tau t-G(\tau)\}.
\end{equation*}
From the relation \eqref{G}, we immediately deduce several inequalities to be utilized later. 
\begin{enumerate}[(a)]
    \item For $t\in [0,\infty)$, $a\in (0,1)$ and $c=p^{1/(q-1)}/q^{1/(p-1)}$,
    \begin{equation}\label{G,a<1}
        \left\{\begin{aligned}
            &a^{q}G(t)\leq G(at)\leq a^{p}G(t), &{} &a^{\frac{1}{p}}G^{-1}(t)\leq G^{-1}(at)\leq a^{\frac{1}{q}}G^{-1}(t),\\
            &\frac{p}{q}a^{q-1}g(t)\leq g(at)\leq \frac{q}{p}a^{p-1}g(t), &{} &ca^{\frac{1}{p-1}}g^{-1}(t)\leq g^{-1}(at)\leq ca^{\frac{1}{q-1}}g^{-1}(t),\\
        \end{aligned}\right.
    \end{equation}
    and for $t\in [0,\infty)$, $a\in(1,\infty)$,
    \begin{equation}\label{G,a>1}
        \left\{\begin{aligned}
        &a^{p}G(t)\leq G(at)\leq a^{q}G(t), &{} &a^{\frac{1}{q}}G^{-1}(t)\leq G^{-1}(at)\leq a^{\frac{1}{p}}G^{-1}(t),\\
            &\frac{q}{p}a^{p-1}g(t)\leq g(at)\leq \frac{p}{q}a^{q-1}g(t), &{} &ca^{\frac{1}{q-1}}g^{-1}(t)\leq g^{-1}(at)\leq ca^{\frac{1}{p-1}}g^{-1}(t).\\
        \end{aligned}\right.
    \end{equation}
     \item The Young inequality with $\epsilon\in (0,1]$:
     \begin{equation}\label{Young}
         t\tau\leq \epsilon^{1-q}G(t)+\epsilon G^{*}(\tau),\;t,\tau\geq0. 
     \end{equation}
     \item For $t,\tau \geq 0$,
     \begin{equation}\label{G*}
         G^{*}(g(t))\leq (q-1)G(t), 
     \end{equation}
     and
     \begin{equation}\label{G triangle}
         2^{-1}(G(t)+G(\tau))\leq G(t+\tau)\leq 2^{q-1}(G(t)+G(\tau)). 
     \end{equation}
\end{enumerate}

Before giving the definition of weak solutions, we shall introduce the notion of Orlicz-Sobolev space. Providing that $N$-function $G$ satisfies the condition \eqref{G}, the Orlicz space $L^{G}(\Omega)$ is described as:
$$L^{G}(\Omega)=\Bigl\{u \text{ is measurable function in }\Omega\colon\int_{\Omega} G(|u(x)|)\differential x<\infty\Bigr\}.$$
The norm of the above space is 
$$\|u\|_{L^{G}(\Omega)}=\inf\Bigl\{\lambda>0\colon\int_{\Omega}G\Bigl(\frac{u(x)}{\lambda}\Bigr)\differential x\leq 1\Bigr\}.$$
We next introduce the fractional Orlicz-Sobolev space $W^{s,G}(\Omega)$
$$W^{s,G}(\Omega)=\Bigl\{u\in L^{G}(\Omega)\colon \int_{\Omega}\int_{\Omega}G\Bigl(\frac{|u(x)-u(y)|}{|x-y|^{s}}\Bigr)\frac{\differential x\differential y}{|x-y|^{N}}<\infty\Bigr\},$$
equipped with the norm
$$\|u\|_{W^{s,G}(\Omega)}=\|u\|_{L^{G}(\Omega)}+[u]_{W^{s,G}(\Omega)},$$
where the semi-norm is defined as
$$[u]_{W^{s,G}(\Omega)}=\inf\Bigl\{\lambda>0\colon\int_{\Omega}\int_{\Omega}G\Bigl(\frac{|u(x)-u(y)|}{\lambda|x-y|^{s}}\Bigr)\frac{\differential x\differential y}{|x-y|^{N}}\leq 1\Bigr\}.$$
For measurable function $u$ in $\realset^{N}$, we define
$$\mathbb{W}^{s,G}(\Omega)=\Bigl\{u|_{\Omega}\in L^{G}(\Omega)\colon\iint_{C_{\Omega}}G\Bigl(\frac{|u(x)-u(y)|}{|x-y|^{s}}\Bigr)\frac{\differential x\differential y}{|x-y|^{N}}<\infty\Bigr\},$$
where $C_{\Omega}\coloneq(\Omega\times \realset^{N})\cup(\realset^{N}\times\Omega)$.

The definition of weak solutions to \cref{luf} is provided as follows:
\begin{definition}
    A measurable function $u\in\mathbb{W}^{s, G}(\Omega)$ is a weak supersolution(subsolution) of \cref{luf}, if 
    \begin{align*}
        &\relphantom{\geq(\leq)} \iint_{C_{\Omega}}g\left(\frac{|u(x)-u(y)|}{|x-y|^{s}}\right)\frac{u(x)-u(y)}{|u(x)-u(y)|}(\varphi(x)-\varphi(y))\frac{K(x,y)}{|x-y|^{s}}\differential x\differential y
        \\ &\geq(\leq)\int_{\Omega}f(x,u(x))\varphi(x)\differential x,\addtocounter{equation}{1}\tag{\theequation}\label{weak}
    \end{align*}
    for all non-negative functions $\varphi\in\mathbb{W}^{s,G}(\Omega)$ with compact support in $\Omega$. $u$ is said to be a weak solution if and only if it is a weak supersolution and a weak subsolution.
\end{definition}
\section{Energy Estimate}\label{section: energy estimate}
We will establish an energy estimate inequality that considers the impact of the inhomogeneous term.
\begin{proposition}\label{proposition: energy estimate}
Assume that $u$ is a weak supersoution to \eqref{luf}, then there exists a constant $\gamma_{*}>0$ depending only on $s,p,q,N,d_{2},\lambda,\Lambda$, such that for any $B_{r}(x_{0})\subset B_{R}(x_{0})\subset\Omega(R<1)$ and any $k\in \realset$, 
\begin{align*}
    &\phantomeq \int_{B_{r}}\int_{B_{r}}G\Bigl(\frac{|\omega_{-}(x)-\omega _{-}(y)|}{|x-y|^{s}}\Bigr)\frac{\differential x\differential y}{|x-y|^{N}}+\int_{B_{r}}\omega_{-}(x)\int_{\mathbb{R}^{N}}g\Bigl(\frac{\omega_{+}(y)}{|x-y|^{s}}\Bigr)\frac{\differential y\differential x}{|x-y|^{N+s}}\\
    &\leq \gamma_{*}\Bigl[G^{*}(d_{1}R^{s})+G\Bigl(\frac{k}{R^{s}}\Bigr)\Bigr]|A_{-}(k,R)|+\gamma_{*}\frac{R^{q}}{(R-r)^{q}}\int_{B_{r}}G\Bigl(\frac{\omega_{-}(x)}{R^{s}}\Bigr)\differential x
    \\&\phantomeq +\gamma_{*}\frac{R^{N+sq}}{(R-r)^{N+sq}}\|\omega_{-}\|_{L^{1}(B_{R})}\Tail(\omega_{-};B_{R}(x_{0})),
\end{align*}
where $\omega_{-}=(u-k)_{-}$ and $A_{-}(k,R)=\{x\in B_{R}\colon u(x)<k\}$.
\end{proposition}
\begin{proof}
    Assume $x_{0}=0$ for simplicity. Let $\omega_{\pm}=(u-k)_{\pm}$. Now we take a cutoff function $\eta\in C^{\infty}_{0}(B_{R})$ with $0\leq\eta\leq1$, vanishing outside $B_{\frac{r+R}{2}}$ and equal to $1$ in $B_{r}$, such that $|\nabla\eta|\leq \frac{4}{R-r}$. Now we select $ \varphi \coloneq\eta^{q}\omega_{-}$ as a test function in the weak formulation \eqref{weak} and obtain that
   \begin{align*}
        &\phantomeq \int_{\realset^{N}}f(x,u(x))\varphi(x)\differential x
        \\&\leq\Lambda\int_{B_{R}}\int_{B_{R}}{g\left(\frac{|u(x)-u(y)|}{|x-y|^{s}}\right)\frac{u(x)-u(y)}{|u(x)-u(y)|}\frac{\varphi(x)-\varphi(y)}{|x-y|^{N+s}}\differential x \differential y}
        \\&\relphantom{=}+2\Lambda\int_{B_{R}}\int_{\realset^{N}\setminus B_{R}}{g\left(\frac{|u(x)-u(y)|}{|x-y|^{s}}\right)\frac{u(x)-u(y)}{|u(x)-u(y)|}\frac{\varphi(x)}{|x-y|^{N+s}}\differential x \differential y}
        \\&\eqcolon\Lambda I_{1}+\Lambda I_{2}.\addtocounter{equation}{1}\tag{\theequation}\label{weak in 3}
    \end{align*}
    
    \noindent\textbf{Estimate of} $\boldsymbol{I_{1}}$. If $x,y\notin A_{-}(k,R)$, then using the fact that $\supp\varphi\subset\subset A_{-}(k,R)$ we know
    \begin{equation}\label{AA}
        g\left(\frac{|u(x)-u(y)|}{|x-y|^{s}}\right)\frac{u(x)-u(y)}{|u(x)-u(y)|}\frac{\varphi(x)-\varphi(y)}{|x-y|^{s}}=0.
    \end{equation}
    
    If $x\in A_{-}(k,R)$, $y\notin A_{-}(k,R)$, then
    \begin{align*}
        &\phantomeq g\left(\frac{|u(x)-u(y)|}{|x-y|^{s}}\right)\frac{u(x)-u(y)}{|u(x)-u(y)|}\frac{\varphi(x)-\varphi(y)}{|x-y|^{s}}
        \\&\leq-\frac{p}{2}\eta^{q}(x)G\left(\frac{\omega_{-}(x)}{|x-y|^{s}}\right)-\frac{1}{2}\eta^{q}(x)\frac{\omega_{-}(x)}{|x-y|^{s}}g\left(\frac{\omega_{+}(y)}{|x-y|^{s}}\right)
        \\&\leq-\frac{p}{2}\min\{\eta^{q}(x),\eta^{q}(y)\}G\left(\frac{|\omega_{-}(x)-\omega_{-}(y)|}{|x-y|^{s}}\right)
        \\&\relphantom{\leq}-\frac{1}{2}\min\{\eta^{q}(x),\eta^{q}(y)\}\frac{\omega_{-}(x)}{|x-y|^{s}}g\left(\frac{\omega_{+}(y)}{|x-y|^{s}}\right).\addtocounter{equation}{1}\tag{\theequation}\label{AB}
    \end{align*}
    
    Inequality \eqref{AB} can also be deduced immediately, since, if $x\in A_{-}(k,R)$, $y\notin A_{-}(k,R)$, 
    \begin{align*}
        g\left(\frac{|u(x)-u(y)|}{|x-y|^{s}}\right)\frac{u(x)-u(y)}{|u(x)-u(y)|}(\varphi(x)-\varphi(y))
        =-\eta^{q}(x)\omega_{-}(x)g\left(\frac{\omega_{-}(x)+\omega_{+}(y)}{|x-y|^{s}}\right).
    \end{align*}
    Observe that
    \begin{equation*}
        g\left(\frac{\omega_{-}(x)+\omega_{+}(y)}{|x-y|^{s}}\right)
        \geq \frac{1}{2}\left[g\left(\frac{\omega_{-}(x)}{|x-y|^{s}}\right)+g\left(\frac{\omega_{+}(y)}{|x-y|^{s}}\right)\right],
    \end{equation*}
    since $g$ is strictly increasing and \eqref{AB} is derived by \eqref{G}.
    
    If $x,y\in A_{-}(k,R)$, then 
    \begin{align*}
        &\phantomeq g\left(\frac{|u(x)-u(y)|}{|x-y|^{s}}\right)\frac{u(x)-u(y)}{|u(x)-u(y)|}\frac{\varphi(x)-\varphi(y)}{|x-y|^{s}}
        \\&\leq -\frac{p}{2}G\left(\frac{|\omega_{-}(x)-\omega_{-}(y)|}{|x-y|^{s}}\right)\min\{\eta^{q}(x),\eta^{q}(y)\}
        \\&\phantomeq +\gamma'G\left(\frac{|\eta(x)-\eta(y)|}{|x-y|^{s}}\max\{\omega_{-}(x),\omega_{-}(y)\}\right),\addtocounter{equation}{1}\tag{\theequation}\label{BB}
    \end{align*}
    where the constant $\gamma'$ depends on $p,q$. 
    Finally, we assume without loss of generality that $k\geq u(x)\geq u(y)$. Considering the fact that $x,y\in A_{-}(k,R)$, we have
    \begin{align*}
        &\phantomeq g\left(\frac{|u(x)-u(y)|}{|x-y|^{s}}\right)\frac{u(x)-u(y)}{|u(x)-u(y)|}\frac{\varphi(x)-\varphi(y)}{|x-y|^{s}}
        \\&=g\left(\frac{|\omega_{-}(x)-\omega_{-}(y)|}{|x-y|^{s}}\right)\frac{\eta^{q}(x)\omega_{-}(x)-\eta^{q}(y)\omega_{-}(y)}{|x-y|^{s}}\eqcolon K.
    \end{align*}
    Notice that, if $\eta(x)\leq\eta(y)$, we have
    \begin{equation*}
        \begin{aligned}
            K&\leq -g\left(\frac{|\omega_{-}(x)-\omega_{-}(y)|}{|x-y|^{s}}\right)\frac{|\omega_{-}(x)-\omega_{-}(y)|}{|x-y|^{s}}\eta^{q}(y)
            \\&\leq-p\eta^{q}(y)G\left(\frac{|\omega_{-}(x)-\omega_{-}(y)|}{|x-y|^{s}}\right)
            \\&\leq -pG\left(\frac{|\omega_{-}(x)-\omega_{-}(y)|}{|x-y|^{s}}\right)\min\{\eta^{q}(x),\eta^{q}(y)\}.
        \end{aligned}
    \end{equation*}
    If $\eta(x)>\eta(y)$, then
    \begin{equation*}
        \begin{aligned}
            K&=-g\left(\frac{|\omega_{-}(x)-\omega_{-}(y)|}{|x-y|^{s}}\right)\frac{|\omega_{-}(x)-\omega_{-}(y)|\eta^{q}(x)}{|x-y|^{s}}
            \\&\relphantom{=}+g\left(\frac{|\omega_{-}(x)-\omega_{-}(y)|}{|x-y|^{s}}\right)\frac{\omega_{-}(y)(\eta^{q}(x)-\eta^{q}(y))}{|x-y|^{s}}\eqcolon K_{1}+K_{2}.
        \end{aligned}
    \end{equation*}
    Using \eqref{G}, we have
    \begin{equation}\label{I1}
        \begin{aligned}
            K_{1}\leq -p\eta^{q}(x)G\left(\frac{|\omega_{-}(x)-\omega_{-}(y)|}{|x-y|^{s}}\right).
        \end{aligned}
    \end{equation}
    Observing the fact that $\eta^{q}(x)-\eta^{q}(y)\leq q \eta^{q-1}(x)(\eta(x)-\eta(y))$, we further compute
    \begin{equation}\label{I2start}
    \begin{aligned}
        K_{2}\leq qg\left(\frac{|\omega_{-}(x)-\omega_{-}(y)|}{|x-y|^{s}}\right)\eta^{q-1}(x)\frac{\eta(x)-\eta(y)}{|x-y|^{s}}\omega_{-}(y).
    \end{aligned}
    \end{equation}
    Then, we apply \eqref{Young} and \eqref{G*} with $\epsilon=\min\{\frac{p}{2(q-1)}, \frac{1}{2}\}$ to obtain that
    \begin{align*}
        &\phantomeq g\left(\frac{|\omega_{-}(x)-\omega_{-}(y)|}{|x-y|^{s}}\right)\eta^{q-1}(x)\frac{\eta(x)-\eta(y)}{|x-y|^{s}}\omega_{-}(y)
        \\&\leq \epsilon G^{*}\left(g\left(\frac{|\omega_{-}(x)-\omega_{-}(y)|}{|x-y|^{s}}\right)\eta^{q-1}(x)\right)+\gamma(\epsilon)G\left(\frac{\eta(x)-\eta(y)}{|x-y|^{s}}\omega_{-}(y)\right)
        \\&\leq\frac{p}{2}G\left(\frac{|\omega_{-}(x)-\omega_{-}(y)|}{|x-y|^{s}}\right)\eta^{q}(x)+\gamma(p,q)G\left(\frac{\eta(x)-\eta(y)}{|x-y|^{s}}\omega_{-}(y)\right).\addtocounter{equation}{1}\tag{\theequation}\label{I2end}
    \end{align*}
    By means of \eqref{I1}--\eqref{I2end}, we finally derive \eqref{BB}.
    Utilizing \eqref{AA}--\eqref{BB} and considering the fact that $\eta=1$ in $B_{r}$, we derive
    \begin{align*}
        &\phantomeq \iint_{B_{R}\times B_{R}}g\left(\frac{|u(x)-u(y)|}{|x-y|^{s}}\right)\frac{u(x)-u(y)}{|u(x)-u(y)|}\frac{\varphi(x)-\varphi(y)}{|x-y|^{N+s}}\differential x\differential y
        \\&\leq \gamma(p,q)\iint_{B_{R}\times B_{R}}G\left(\frac{|\eta(x)-\eta(y)|}{|x-y|^{s}}\max\{\omega_{-}(x),\omega_{-}(y)\}\right)\frac{1}{|x-y|^{N}}\differential x\differential y
        \\&\phantomeq -\frac{p}{2}\iint_{B_{r}\times B_{r}}G\left(\frac{|\omega_{-}(x)-\omega_{-}(y)|}{|x-y|^{s}}\right)\frac{1}{|x-y|^{N}}\differential x\differential y
        \\&\phantomeq -\frac{1}{2}\iint_{B_{r}\times B_{r}}\frac{\omega_{-}(x)}{|x-y|^{N+s}}g\left(\frac{\omega_{+}(y)}{|x-y|^{s}}\right)\differential x \differential y.\addtocounter{equation}{1}\tag{\theequation}\label{BRBR}
    \end{align*}
    Recalling that $|\nabla\eta|\leq \frac{4}{R-r}$, we may apply \eqref{G,a>1} and Mean Value Theorem to yield 
    \begin{align*}
        &\phantomeq \int_{B_{R}}\int_{B_{R}}G\left(\frac{|\eta(x)-\eta(y)|}{|x-y|^{s}}\max\{\omega_{-}(x),\omega_{-}(y)\}\right)\frac{1}{|x-y|^{N}}\differential x\differential y
        \\&\leq2\int_{B_{R}}\int_{B_{R}}G\left(\frac{|\eta(x)-\eta(y)|}{|x-y|^{s}}\omega_{-}(x)\right)\frac{1}{|x-y|^{N}}\differential x\differential y
        \\&\leq2\int_{B_{R}}\int_{B_{R}}G\left(\frac{|\nabla\eta||x-y|}{|x-y|^{s}}\omega_{-}(x)\right)\frac{1}{|x-y|^{N}}\differential x\differential y
        \\&\leq 2\int_{B_{R}}\int_{B_{R}}G\left(\frac{8R}{R-r}\frac{(2R)^{s}}{|x-y|^{s}}\frac{\omega_{-}(x)}{(2R)^{s}}\right)\frac{1}{|x-y|^{N}}\differential x\differential y
        \\&\leq \gamma(N,s,p,q)\frac{R^{sq+q}}{(R-r)^{q}}\int_{B_{R}}G\left(\frac{\omega_{-}(x)}{R^{s}}\right)\int_{B_{R}}\frac{1}{|x-y|^{N+sq}}\differential y\differential x
        \\&\leq \gamma(N,s,p,q) \frac{R^{q}}{(R-r)^{q}}\int_{B_{R}}G\left(\frac{\omega_{-}(x)}{R^{s}}\right)\differential x,
    \end{align*}
    in which the penultimate inequality arises from the facts $\frac{R}{R-r}>1$ and $\frac{2R}{|x-y|}>1$. Substituting the above inequality to \eqref{BRBR} shows that 
    \begin{align*}
        &\phantomeq \int_{B_{R}}\int_{B_{R}}g\left(\frac{|u(x)-u(y)|}{|x-y|^{s}}\right)\frac{u(x)-u(y)}{|u(x)-u(y)|}\frac{\varphi(x)-\varphi(y)}{|x-y|^{N+s}}\differential x\differential y
        \\&\leq-\frac{p}{2}\int_{B_{r}}\int_{B_{r}}G\left(\frac{|\omega_{-}(x)-\omega_{-}(y)|}{|x-y|^{s}}\right)\frac{1}{|x-y|^{N}}\differential x\differential y
        \\&\phantomeq -\frac{1}{2}\int_{B_{r}}\int_{B_{r}}\frac{\omega_{-}(x)}{|x-y|^{N+s}}g\left(\frac{\omega_{+}(y)}{|x-y|^{s}}\right)\differential x \differential y
        \\&\phantomeq +\gamma(N,s,p,q)\frac{R^{q}}{(R-r)^{q}}\int_{B_{R}}G\left(\frac{\omega_{-}(x)}{R^{s}}\right)\differential x.\addtocounter{equation}{1}\tag{\theequation}\label{BRBR2}
    \end{align*}
    Next, we will deal with the second term on the right-hand side of \eqref{weak in 3}.
    
    \noindent\textbf{Estimate of} $\boldsymbol{I_{2}}$.    
    We now turn our attention to the term integrated on $B_{R}\times(\realset^{N}\setminus B_{R})$.
    \begin{equation*}
        \begin{aligned}
            &\phantomeq 2\int_{B_{R}}\eta^{q}(x)\omega_{-}(x)\int_{\realset^{N}\setminus B_{R}}g\left(\frac{|u(x)-u(y)|}{|x-y|^{s}}\right)\frac{u(x)-u(y)}{|u(x)-u(y)|}\frac{1}{|x-y|^{N+s}}\differential y\differential x
            \\&\leq 2\int_{A_{-}(k,\frac{r+R}{2})}\omega_{-}(x)\int_{\realset^{N}\cap \{u(y)<u(x)\}\setminus B_{R}}g\left(\frac{|u(x)-u(y)|}{|x-y|^{s}}\right)\frac{1}{|x-y|^{N+s}}\differential y\differential x
            \\&\relphantom{\leq}-2\int_{A_{-}(k,r)}\omega_{-}(x)\int_{\realset^{N}\cap\{u(y)\geq u(x)\}\setminus B_{R}}g\left(\frac{|u(x)-u(y)|}{|x-y|^{s}}\right)\frac{1}{|x-y|^{N+s}}\differential y\differential x
            \\&\eqcolon J_{1}-J_{2}.
        \end{aligned}
    \end{equation*}
    Given $x\in B_{\frac{r+R}{2}}$ and $y\in \realset^{N}\setminus B_{R}$, we have $$|x-y|\geq \frac{R-r}{2R}|y|.$$
    This together with \eqref{G,a>1} result in
    \begin{align*}
        J_{1}&\leq 2\int_{B_{\frac{r+R}{2}}}\omega_{-}(x)\int_{\realset^{N}\setminus B_{R}}g\left(\frac{\omega_{-}(y)}{|y|^{s}}\left(\frac{2R}{R-r}\right)^{s}\right)\frac{1}{|y|^{N+s}}\left(\frac{2R}{R-r}\right)^{N+s}\differential y\differential x
        \\&\leq \gamma(N,s,p,q)\frac{R^{N+sq}}{(R-r)^{N+sq}}
        \int_{B_{\frac{r+R}{2}}}\omega_{-}(x)\int_{\realset^{N}\setminus B_{R}}g\left(\frac{\omega_{-}(y)}{|y|^{s}}\right)\frac{1}{|y|^{N+s}}\differential y\differential x
        \\&\leq \gamma(N,s,p,q)\frac{R^{N+sq}}{(R-r)^{N+sq}}\|\omega_{-}(x)\|_{L^{1}(B_{R})}\Tail(\omega_{-};B_{R}).\addtocounter{equation}{1}\tag{\theequation}\label{J1}
    \end{align*}
    It's easy to check that
    \begin{equation*}
        \{u(y)\geq k\}\subset \{u(y)\geq u(x)\},
    \end{equation*}
    which leads to
    \begin{align*}
        J_{2}&\geq 2\int_{B_{r}}\omega_{-}(x)\int_{\realset^{N}\cap A_{+}(k)\setminus B_{R}}g\left(\frac{\omega_{+}(y)+\omega_{-}(x)}{|x-y|^{s}}\right)\frac{1}{|x-y|^{N+s}}\differential y\differential x
        \\&\geq 2\int_{B_{r}}\omega_{-}(x)\int_{\realset^{N}\cap A_{+}(k)\setminus B_{R}}g\left(\frac{\omega_{+}(y)}{|x-y|^{s}}\right)\frac{1}{|x-y|^{N+s}}\differential y\differential x
        \\&=2\int_{B_{r}}\omega_{-}(x)\int_{\realset^{N}\setminus B_{R}}g\left(\frac{\omega_{+}(y)}{|x-y|^{s}}\right)\frac{1}{|x-y|^{N+s}}\differential y\differential x,\addtocounter{equation}{1}\tag{\theequation}\label{J2}
    \end{align*}
    where $A_{+}(k)=\{x\in \realset^{N} \colon u(x)<k\}$. 
    Here, to obtain the last line we use the fact that $\omega_{+}(y)=0,$ while $y\notin \realset^{N}\cap A_{+}(k)$. It follows from \eqref{J1} and \eqref{J2} that
    \begin{align*}
        % \renewcommand{\indent}{\expandafter\expandafter\iint}
        &\phantomeq \int_{B_{R}}\int_{\realset^{N}\setminus B_{R}}g\left(\frac{|u(x)-u(y)|}{|x-y|^{s}}\right)\frac{u(x)-u(y)}{|u(x)-u(y)|}\frac{\varphi(x)-\varphi(y)}{|x-y|^{N+s}}\differential x\differential y
        \\&\leq \gamma_{1}(N,s,p,q)\frac{R^{N+sq}}{(R-r)^{N+sq}}\|\omega_{-}(x)\|_{L^{1}(B_{R})}\Tail(\omega_{-};B_{R})
        \\&\phantomeq -\gamma_{2}(N,s,p,q)\int_{B_{r}}\omega_{-}(x)\int_{\realset^{N}\setminus B_{R}}g\left(\frac{\omega_{+}(y)}{|x-y|^{s}}\right)\frac{1}{|x-y|^{N+s}}\differential y\differential x.\addtocounter{equation}{1}\tag{\theequation}\label{BRRnBR}
    \end{align*}
    To finish the proof, we shall estimate the term on the left-hand side of \eqref{weak in 3}. Recalling the property of $f(x,u(x))$, we get
    \begin{equation}\label{f}
        \int_{\realset^{N}}|f(x,u(x))|\varphi(x)\differential x\leq \int_{A_{-}(k,\frac{r+R}{2})}\bigl(d_{1}+d_{2}g(u(x))\bigr)\eta^{q}(x)\omega_{-}(x)\differential x.
    \end{equation}
    For the first term of the inequality above, we may apply the Young inequality \eqref{Young} with $\epsilon=\left(\frac{R}{R-r}\right)^{\frac{q}{1-q}}<1$, to obtain
    \begin{align*}
        &\phantomeq \int_{A_{-}(k,\frac{r+R}{2})}d_{1}\eta^{q}(x)\omega_{-}(x)\differential x
        \\&\leq \left(\frac{R}{R-r}\right)^{\frac{q}{1-q}}G^{*}(d_{1}R^{s})|A_{-}(k,R)|+\left(\frac{R}{R-r}\right)^{q}\int_{B_{R}}G\left(\frac{\omega_{-}(x)}{R^{s}}\right)\differential x.\addtocounter{equation}{1}\tag{\theequation}\label{f1}
    \end{align*}
    For the second term of the inequality \eqref{f}, we notice that if $x\in A_{-}(k)$
    \begin{align*}
        g(|u(x)|)=g(|\omega_{-}(x)-k|)\leq g(\omega_{-}(x)+k)\leq C(p,q) g(\omega_{-}(x))+\gamma(p,q) g(k),
    \end{align*}
    then
    \begin{align*}
        &\phantomeq \int_{A_{-}(k,\frac{r+R}{2})}d_{2}g(|u(x)|)\eta^{q}(x)\omega_{-}(x)\differential x
        \\&\leq \gamma\int_{A_{-}(k,\frac{r+R}{2})}\eta^{q}d_{2}\left(g(\omega_{-}(x)\omega_{-}(x)+g(k)\omega_{-}(x)\right)\differential x
        \\&\leq \widetilde{\gamma}\left(\frac{R}{R-r}\right)^{q}\int_{B_{R}}G\left(\frac{\omega_{-}(x)}{R^{s}}\right)\differential x + \widetilde{\gamma}G\left(\frac{k}{R^{s}}\right)|A_{-}(k,R)|,\addtocounter{equation}{1}\tag{\theequation}\label{f2}
    \end{align*}
    where $\widetilde{\gamma}$ depends on the data $\{N,s,p,q,d_{2}\}$ and we use the fact $R<1$. 
    By putting together \eqref{BRBR2}, \eqref{BRRnBR}, \eqref{f1} and \eqref{f2} and adjusting the constants, we finally complete the proof.
\end{proof}
\section{DENSITY LEMMA}\label{section: density lemma}

The following lemma shows the spread of pointwise positivity in space. 
 \begin{lemma}\label{lemma1}
        Let $k\geq 0$ and $R\leq 1$ be parameters. Assume that $u$ is a supersolution of \eqref{luf}, non-negative in $B_R\left(x_0\right)\subseteq \Omega$. Then there exists a constant $\nu\in (0,1)$ depending only on the data $\left\{s,p,q,N,d_{1},d_{2}\right\}$, such that if
        \begin{equation*}
            \left| \left\{u\leq k\right\} \cap B_r\left(x_0\right) \right| \leq \nu \left|B_r\left(x_0\right)\right|,
        \end{equation*}
        then either 
        \begin{equation*}
            r^s g^{-1}\left(r^s \Tail\left(u_{-};B_R\left(x_0\right)\right)\right)+r^sG^{-1}G^{*}\left(d_1 r^s\right)>k,
        \end{equation*}
        or
        \begin{equation*}
            u\geq \frac{1}{2}k, \forall x\in B_{\frac{1}{2}r}\left(x_0\right),
        \end{equation*}
        where $B_r\left(x_0\right)\subseteq B_R\left(x_0\right)$.
    \end{lemma}
    \begin{proof}
    Without loss of generality, we suppose $x_{0}=0$ and 
    \begin{equation}\label{assumption}
        r^{s}g^{-1}\left(r^{s}\Tail(u_{-};B_{R})\right)+r^{s}G^{-1}G^{*}(d_{1}r^{s})\leq k.
    \end{equation}
    For all $n\in \mathbb{N}_{0}$, set
    \begin{equation*}
    \left\{
    \begin{gathered}
        k_{n}=\frac{k}{2}+\frac{k}{2^{n+1}},\\
        r_{n}=\frac{r}{2}+\frac{r}{2^{n+1}},\,\tilde{r}_{n}=\frac{r_{n}+r_{n+1}}{2},\\
        \hat{r}_{n}=\frac{3r_{n}+r_{n+1}}{4},\,\overline{r}_{n}=\frac{r_{n}+3r_{n+1}}{4},\\
        B_{n}=B_{r_{n}},\,\tilde{B}_{n}=B_{\tilde{r}_{n}},\,\hat{B}_{n}=B_{\hat{r}_{n}},\,\overline{B}_{n}=B_{\overline{r}_{n}}.
    \end{gathered}
    \right.
    \end{equation*}
    Observe that $B_{n+1}\subset\overline{B}_{n}\subset \tilde{B}_{n}\subset\hat{B}_{n}\subset B_{n}$. Now we take a cutoff function $\phi$ in $B_{n}$, vanishing outside $\hat{B}_{n}$, and equal to the identity in $\tilde{B}_{n}$, such that $$|D\phi|\leq \frac{2^{n}}{r}.$$
   Selecting $k=k_{n},B_{r}=\tilde{B}_{n}\:\text{and}\:B_{R}=B_{n}$, we can make use of the energy estimate of Proposition \ref{proposition: energy estimate}. As a result, we have
    \begin{align*}
        &\phantomeq \int_{B_{n+1}}\int_{B_{n+1}}G\left(\frac{|\omega_{-}(x)-\omega_{-}(y)|}{|x-y|^{s}}\right)\frac{\differential x\differential y}{|x-y|^{N}}
        \\&\leq \gamma_{*}\Bigl[G^{*}(d_{1}r^{s})+G\Bigl(\frac{k_{n}}{r^{s}}\Bigr)\Bigr]|B_{n}\cap\{u<k_{n}\}|+\gamma_{*}\frac{r_{n}^{q}}{(r_{n}-\tilde{r}_{n})^{q}}\int_{B_{n}}G\Bigl(\frac{\omega_{-}(x)}{r_{n}^{s}}\Bigr)\differential x
        \\&\phantomeq +\gamma_{*}\frac{r_{n}^{N+sq}}{(r_{n}-\tilde{r}_{n})^{N+sq}}\|\omega_{-}\|_{L^{1}(B_{n})}\Tail(\omega_{-};B_{n}),\addtocounter{equation}{1}\tag{\theequation}\label{equation: estimate}
    \end{align*}
    where $\omega_{-}=(u-k_{n})_{-}$. Firstly we shall focus on estimating the first term on the right-hand side of \eqref{equation: estimate}. Recalling \eqref{assumption} that $r^{s}G^{-1}G^{*}(d_{1}r^{s})\leq k$, we arrive that
    \begin{equation}\label{equation: the 1st term}
        \begin{aligned}
            \left[G^{*}(d_{1}r^{s})+G\left(\frac{k_{n}}{r^{s}}\right)\right]|B_{n}\cap \{u<k_{n}\}|\leq 2G\left(\frac{k_{n}}{r^{s}}\right)|B_{n}\cap \{u<k_{n}\}|.
        \end{aligned}
    \end{equation}
    Considering the definitions of $r_{n}$ and $\tilde{r}_{n}$, there holds 
    \begin{equation*}
        \frac{r_{n}}{r_{n}-\tilde{r}_{n}}=2^{n+2}+4.
    \end{equation*}
    Consequently, we estimate the second term as
    \begin{align*}
        \frac{r_{n}^{q}}{(r_{n}-\tilde{r}_{n})^{q}}\int_{B_{n}}G\Bigl(\frac{\omega_{-}(x)}{r_{n}^{s}}\Bigr)\differential x
        &\leq \gamma(p,q)2^{nq}\int_{B_{n}}G\left(\frac{\omega_{-}(x)}{r^{s}}\right)\differential x
        \\&\leq \gamma(p,q)2^{nq}G\left(\frac{k_{n}}{r^{s}}\right)|B_{n}\cap\{u<k_{n}\}|.\addtocounter{equation}{1}\tag{\theequation}\label{equation: the 2nd term}
    \end{align*}
    The last term in \eqref{equation: estimate} shall be estimated as follows. An application of \eqref{G triangle} and \eqref{G} gives that
    \begin{equation*}
        g\left(\frac{k_{n}+u_{-}}{|x_{0}-y|^{s}}\right)\leq \frac{q2^{q-1}}{p}\left[g\left(\frac{k_{n}}{|x_{0}-y|^{s}}\right)+g\left(\frac{u_{-}(y)}{|x_{0}-y|^{s}}\right)\right].
    \end{equation*}
    We further get
    \begin{align*}
        &\phantomeq \frac{r_{n}^{N+sq}}{(r_{n}-\tilde{r}_{n})^{N+sq}}\|\omega_{-}\|_{L^{1}(B_{n})}\Tail(\omega_{-};B_{n})
        \\&\leq \gamma(p,q)2^{n(N+sq)}k_{n}|B_{n}\cap\{u<k_{n}\}|\int_{\realset^{N}\setminus B_{n}}g\left(\frac{k_{n}+u_{-}(y)}{|x_{0}-y|^{s}}\right)\frac{\differential y}{|x_{0}-y|^{N+s}}
        \\&\leq \gamma(p,q)2^{n(N+sq)}k_{n}|B_{n}\cap \{u<k_{n}\}|\int_{\realset^{N}\setminus B_{n}}\left[g\left(\frac{k_{n}}{r_{n}^{s}}\right)+g\left(\frac{u_{-}(y)}{|x_{0}-y|^{s}}\right)\right]\frac{\differential y}{|x_{0}-y|^{N+s}}
        \\&\leq \gamma(p,q)2^{n(N+sq)}k_{n}|B_{n}\cap \{u<k_{n}\}|\left[r^{-s}g\left(\frac{k_{n}}{r^{s}}\right)+\Tail(u_{-};B_{R})\right].
    \end{align*}
    By means of \eqref{G} and \eqref{assumption} $r^{s}g^{-1}(r^{s}\Tail(u_{-};B_{R}))\leq k$, we arrive at
    \begin{equation}\label{equation: the 3rd term}
        \frac{r_{n}^{N+sq}}{(r_{n}-\tilde{r}_{n})^{N+sq}}\|\omega_{-}\|_{L^{1}(B_{n})}\Tail(\omega_{-};B_{n})\leq \gamma2^{n(N+sq)}k_{n}\ |B_{n}\cap \{u<k_{n}\}|\ r^{-s}g\Bigl(\frac{k_{n}}{r^{s}}\Bigr).
    \end{equation}
    Combining \eqref{equation: estimate}--\eqref{equation: the 3rd term}, we can get
    \begin{align*}
        &\phantomeq \int_{B_{n+1}}\int_{B_{n+1}}G\left(\frac{|\omega_{-}(x)-\omega_{-}(y)|}{|x-y|^{s}}\right)\frac{\differential x\differential y}{|x-y|^{N}}
        \\&\leq \gamma(p,q) 2^{n(N+sq+q)}G\left(\frac{k_{n}}{r^{s}}\right)|B_{n}\cap \{u<k_{n}\}|.
    \end{align*}
    According to Lemma 4.1 in \cite{byun2021local}, there exists a constant $\theta>1$ depending only on $N,s$, such that
    \begin{align*}
        &\phantomeq \left[\bbint{B_{n+1}}G^{\theta}\left(\frac{|\omega_{-}-(\omega_{-})_{B_{n+1}}|}{r_{n+1}^{s}}\right)\differential x\right]^{\frac{1}{\theta}}
        \\&\leq \gamma(p,q)\bbint{B_{n+1}}\int_{B_{n+1}}G\left(\frac{|\omega_{-}(x)-\omega_{-}(y)|}{|x-y|^{s}}\right)\frac{\differential x\differential y}{|x-y|^{N}}\notag
        \\&\leq \gamma(p,q)2^{n(N+sq+q)}G\left(\frac{k_{n}}{r^{s}}\right)\frac{|B_{n}\cap \{u<k_{n}\}|}{|B_{n}|}.
    \end{align*}
    We further employ the Jensen inequality to get the following display
    \begin{align*}
        &\phantomeq \left[\bbint{B_{n+1}}G^{\theta}\left(\frac{\omega_{-}}{r^{s}_{n+1}}\right)\differential x\right]^{\frac{1}{\theta}}
        \\ &\leq \gamma\bbint{B_{n+1}}G\left(\frac{\omega_{-}}{r_{n+1}^{s}}\right)\differential x
        +\gamma\left[\bbint{B_{n+1}}G^{\theta}\left(\frac{|\omega_{-}-(\omega_{-})_{B_{n+1}}|}{r_{n+1}^{s}}\right)\differential x\right]^{\frac{1}{\theta}}.
    \end{align*}
    Recall the definition of $r_{n+1}$ with $\frac{r}{2}\leq r_{n+1}\leq r, $ then
    $$G^{\theta}\left(\frac{\omega_{-}}{r_{n+1}^{s}}\right)\geq G^{\theta}\left(\frac{k_{n}-k_{n+1}}{r^{s}_{n+1}}\right)\chi_{\{u<k_{n+1}\}}\geq \gamma'2^{-nq\theta}G\left(\frac{k_{n}}{r^{s}}\right)\chi_{\{u<k_{n+1}\}}.$$
    It can be deduced that
    \begin{equation*}
        \gamma2^{-jq}G\left(\frac{k_{n}}{r^{s}}\right)\left(\bbint{B_{j+1}}\chi_{\{u<k_{j+1} \}}\differential x\right)^{\frac{1}{\theta}}\leq
        \gamma 2^{n(N+sq+q)}G\left(\frac{k_{n}}{r^{s}}\right)\frac{|B_{n}\cap \{u<k_{n}\}|}{|B_{n}|},
    \end{equation*}
    which means
    \begin{equation*}
        \begin{aligned}
            \left(\frac{|B_{n+1}\cap\{u<k_{n+1}\}|}{|B_{n+1}|}\right)^{\frac{1}{\theta}}\leq \gamma 2^{n(N+sq+2q)}\frac{|B_{n}\cap\{u<k_{n}\}|}{|B_{n}|}.
        \end{aligned}
    \end{equation*}
    Denote $$A_{n}\coloneq\frac{|B_{n}\cap \{u<k_{n}\}|}{B_{n}}.$$ Then we have
    $$A_{n+1}\leq \gamma 2^{n(N+sq+2q)\theta}A_{n}^{\theta}.$$
    According to the Lemma 7.1 in \cite{giusti2003direct}, we have to make sure that
    \begin{align*}
        A_{0}=\frac{|B_{r}\cap \{u<k\}|}{|B_{r}|}\leq \gamma^{\frac{-1}{\theta-1}}2^{-(N+sq+2q)\frac{\theta}{(\theta-1)^{2}}}
    \end{align*}
    to obtain $A_{n}\rightarrow 0$ as $n\rightarrow \infty$. Notice that
    $\gamma^{\frac{-1}{\theta-1}}2^{-(N+sq+2q)\frac{\theta}{(\theta-1)^{2}}}<1$ depending only on $\{N,s,p,q,d_{2},\lambda,\Lambda\}$, so there exists a constant $\nu\in(0,1)$ depending on $\{N,s,p,q,d_{2},\lambda,\Lambda\}$ to justify the desired result
    \begin{equation*}
        \lim_{n\rightarrow \infty }A_{n}=0.
    \end{equation*}
    In other words, we draw a conclusion
    \begin{equation*}
        u\geq \frac{k}{2}
    \end{equation*}
    in $B_{\frac{r}{2}}$. We now complete the proof. 
\end{proof}

\section{PROOF OF MAIN THEOREM}
The following two lemmas apply the energy estimate to obtain two different measure density inequalities, which play a key role in proving our main results Theorem \ref{theorem: 2} and Remark \ref{theorem: 1}. 
The first lemma measures the shrinking set with a nonlocal integral.
\begin{lemma}\label{lemma3}
        Let $k\geq 0$ and $R\leq 1$. Let $u\in \mathbb{W}^{s,G}(\Omega)$ be a non-negative supersolution of \eqref{luf} in $B_R\left(x_0\right)\subseteq \Omega$. There exists a constant $\gamma_2>1$ depending only on $\left\{s,p,q,N,d_2,\lambda, \Lambda\right\}$, such that either
        \begin{equation*}
            r^s g^{-1}\left(r^s \Tail\left(u_{-};B_R\left(x_0\right)\right)\right)+r^sG^{-1}G^{*}\left(d_1 r^s\right)>k,
        \end{equation*}
        or
        \begin{equation*}
            \left| \left\{u\leq k\right\} \cap B_r\left(x_0\right) \right| \leq \frac{\gamma_2\,g\left(\frac{k}{r^s}\right)}{r^s \Tail(u_+;B_r\left(x_0\right))} \left|B_r\right|,
        \end{equation*}
        where $B_{2r}\left(x_0\right)\subseteq B_R\left(x_0\right)$.
    \end{lemma}
    \begin{proof}
Without sacrificing generality, we also assume that $x_{0}=0$, and
    $$r^{s}g^{-1}(r^{s}\Tail(u_{-};B_{R}))+r^{s}G^{-1}G^{*}(d_{1}r^{s})\leq k.$$
    Using the similar proceed as in Lemma \ref{lemma2}, we can get
    \begin{equation}\label{equation: lemma estimate 3}
        \int_{B_{r}}\omega_{-}(x)\int_{\realset^{N}}g\left(\frac{\omega_{+}(y)}{|x-y|^{s}}\right)\frac{1}{|x-y|^{N+s}}\differential y\differential x\leq \gamma \frac{k}{r^{s}}g\left(\frac{k}{r^{s}}\right)|B_{r}|,
    \end{equation}
    and
    \begin{equation*}
        g\left(\frac{u_{+}(y)}{|x-y|^{s}}\right)\leq g\left(\frac{(u(y)-k)_{+}+k}{|x-y|^{s}}\right)\leq C\left[g\left(\frac{\omega_{+}(y)}{|x-y|^{s}}\right)+g\left(\frac{k}{|x-y|^{s}}\right)\right].
    \end{equation*}
   Then, we know
    \begin{align*}
        &\phantomeq \int_{B_{r}}\omega_{-}(x)\int_{\realset^{N}}g\left(\frac{\omega_{+}(y)}{|x-y|^{s}}\right)\frac{1}{|x-y|^{N+s}}\differential y\differential x
        \\&\geq \int_{B_{r}}\int_{\realset^{N}\setminus B_{r}}\omega_{-}(x)g\left(\frac{\omega_{+}(y)}{|x-y|^{s}}\right)\frac{1}{|x-y|^{N+s}}\differential y\differential x
        \\&\geq C\int_{B_{r}}\int_{\realset^{N}\setminus B_{r}}\frac{\omega_{-}(x)}{(2|y|)^{N+s}}\left[g\left(\frac{u_{+}(y)}{(2|y|)^{s}}\right)-g\left(\frac{k}{(2|y|)^{s}}\right)\right]\differential y\differential x
        \\&\geq C\int_{B_{r}}\int_{\realset^{N}\setminus B_{r}}\frac{\omega_{-}(x)}{|y|^{N+s}}g\left(\frac{u_{+}(y)}{|y|^{s}}\right)\differential y\differential x-\gamma\frac{k}{r^{s}}g\left(\frac{k}{r^{s}}\right)|B_{r}|.
    \end{align*}
    Combine the above inequality and \eqref{equation: lemma estimate 3} to get
    \begin{equation*}
        \int_{B_{r}}\int_{\realset^{N}\setminus B_{r}}\frac{\omega_{-}(x)}{|y|^{N+s}}g\left(\frac{u_{+}(y)}{|y|^{s}}\right)\differential y\differential x\leq \gamma\frac{k}{r^{s}}g\left(\frac{k}{r^{s}}\right)|B_{r}|.
    \end{equation*}
    Besides we can obtain the following result by the method of Lemma \ref{lemma3}:
    \begin{equation*}
        |\{u<\frac{1}{2}k\}\cap B_{r}|\leq \frac{\gamma g\left(\frac{k}{r^{s}}\right)}{r^{s}\Tail(u_{+};B_{r})}|B_{r}|.
    \end{equation*}
    Then the proof is complete.
\end{proof}

The measure shrinking inequality depends on a local integral associated with $g$ in the second lemma.
\begin{lemma}\label{lemma2}
    Let $k\geq 0$ and $R\leq 1$. Let $u\in \mathbb{W}^{s,G}(\Omega)$ be a non-negative supersolution of \eqref{luf} in $B_R\left(x_0\right)\subseteq \Omega$. There exists a constant $\gamma_1>1$ depending only on $\left\{s,p,q,N,d_2,\lambda, \Lambda\right\}$, such that either
        \begin{equation*}
            r^s g^{-1}\left(r^s \Tail\left(u_{-};B_R\left(x_0\right)\right)\right)+r^sG^{-1}G^{*}\left(d_1 r^s\right)>k,
        \end{equation*}
        or
        \begin{equation*}
            \left| \left\{u\leq k\right\} \cap B_r\left(x_0\right) \right| \leq \frac{\gamma_1\,g\left(\frac{k}{r^s}\right)}{\bbint{B_r}g\bigl(\frac{u(x)}{r^s}\bigr)\differential x} \left|B_r\right|,
        \end{equation*}
        where $B_{2r}\left(x_0\right)\subseteq B_R\left(x_0\right)$.
\end{lemma}
\begin{proof}
    Without loss of generality, we suppose $x_{0}=0$, and
    $$r^{s}g^{-1}(r^{s}\Tail(u_{-};B_{R}))+r^{s}G^{-1}G^{*}(d_{1}r^{s})\leq k.$$
    Apply the energy estimate Proposition \ref{proposition: energy estimate} in $B_{r}\subset B_{2r}$ to get
    \begin{align*}
        &\phantomeq \int_{B_{r}}\omega_{-}(x)\int_{\realset^{N}}g\left(\frac{\omega_{+}(y)}{|x-y|^{s}}\right)\frac{1}{|x-y|^{N+s}}\differential y\differential x
        \\&\leq \gamma_{*}\int_{B_{2r}}G\left(\frac{\omega_{-}(x)}{r^{s}}\right)\differential x+\gamma_{*}\|\omega_{-}\|_{L^{1}(B_{2r})}\Tail(\omega_{-};B_{2r})
        \\&\phantomeq +\gamma_{*}\left[G^{*}(d_{1}r^{s})+G\left(\frac{k}{r^{s}}\right)\right]|A_{-}(k,2r)|
        \\&\leq \gamma G\left(\frac{k}{r^{s}}\right)|B_{r}|+\gamma k|B_{r}|\Tail(\omega_{-};B_{2r})+\gamma G^{*}(d_{1}r^{s})|B_{r}|.\addtocounter{equation}{1}\tag{\theequation}\label{equation: lemma estimate}
    \end{align*}
    According to the assumption $r^{s}g^{-1}(r^{s}\Tail(u_{-};B_{R}))\leq k$, we can estimate the tail term as follows
    \begin{align*}
        \Tail(\omega_{-};B_{2r})&=\int_{\realset^{N}\setminus B_{2r}}g\left(\frac{\omega_{-}(x)}{|x|^{s}}\right)\frac{1}{|x|^{N+s}}\differential x
        \\&=\int_{B_{R}\setminus B_{2r}}g\left(\frac{\omega_{-}(x)}{|x|^{s}}\right)\frac{1}{|x|^{N+s}}\differential x+\int_{\realset^{N}\setminus B_{R}}g\left(\frac{\omega_{-}(x)}{|x|^{s}}\right)\frac{1}{|x|^{N+s}}\differential x
        \\&\leq \gamma r^{-s}g\left(\frac{k}{r^{s}}\right).
    \end{align*}
    Due to the estimate $r^{s}G^{-1}G^{*}(d_{1}r^{s})\leq k$, \eqref{equation: lemma estimate} turns into
    \begin{equation}\label{equation: lemma estimate 2}
        \int_{B_{r}}\omega_{-}(x)\int_{\realset^{N}}g\left(\frac{\omega_{+}(y)}{|x-y|^{s}}\right)\frac{1}{|x-y|^{N+s}}\differential y\differential x\leq \gamma \frac{k}{r^{s}}g\left(\frac{k}{r^{s}}\right)|B_{r}|.
    \end{equation}
    In the following, we will focus on the left part of the above inequality. From \eqref{G} and \eqref{G triangle}, we know
    \begin{equation*}
        g\left(\frac{u_{+}(y)}{|x-y|^{s}}\right)\leq g\left(\frac{(u(y)-k)_{+}+k}{|x-y|^{s}}\right)\leq 2\left[g\left(\frac{\omega_{+}(y)}{|x-y|^{s}}\right)+g\left(\frac{k}{|x-y|^{s}}\right)\right].
    \end{equation*}
    Based on the above inequality, we can see
    \begin{align*}
        &\phantomeq \int_{B_{r}}\omega_{-}(x)\int_{\realset^{N}}g\left(\frac{\omega_{+}(y)}{|x-y|^{s}}\right)\frac{1}{|x-y|^{N+s}}\differential y\differential x
        \\&\geq \int_{B_{r}}\int_{B_{r}}\omega_{-}(x)g\left(\frac{\omega_{+}(y)}{|x-y|^{s}}\right)\frac{1}{|x-y|^{N+s}}\differential y\differential x
        \\&\geq C\int_{B_{r}}\int_{B_{r}}\frac{\omega_{-}(x)}{|x-y|^{N+s}}\left[g\left(\frac{u(y)}{|x-y|^{s}}\right)-g\left(\frac{k}{|x-y|^{s}}\right)\right]\differential y\differential x
        \\&\geq Cr^{-(N+s)}\int_{B_{r}}\int_{B_{r}}\omega_{-}(x)g\left(\frac{u(y)}{r^{s}}\right)\differential y\differential x-\gamma\frac{k}{r^{s}}g\left(\frac{k}{r^{s}}\right)|B_{r}|.
    \end{align*}
    Combine this with \eqref{equation: lemma estimate 2} to reach
    \begin{equation*}
        \int_{B_{r}}\omega_{-}(x)\bbint{B_{r}}g\left(\frac{u(y)}{r^{s}}\right)\differential y\differential x\leq \gamma kg\left(\frac{k}{r^{s}}\right)|B_{r}|.
    \end{equation*}
    Noticing that
    \begin{equation*}
        \int_{B_{r}}\omega_{-}(x)\differential x\geq \frac{1}{2}k|\{u<\frac{1}{2}k\}\cap B_{r}|,
    \end{equation*}
    we arrive that
    \begin{equation*}
        |\{u<\frac{1}{2}k\}\cap B_{r}|\leq \frac{\gamma g\left(\frac{k}{r^{s}}\right)}{\bbint{_{B_{r}}}g\left(\frac{u(x)}{r^{s}}\right)\differential x}|B_{r}|.
    \end{equation*}
    In this way, we complete our proof.
\end{proof}

Next, we present the proof of \cref{theorem: 2} and Remark \ref{theorem: 1}.
\begin{proof}[\textbf{Proof of \ref{theorem: 2} and \ref{theorem: 1}}]
    For simplicity, we choose to omit the symbol $x_{0}$. From Lemma \ref{lemma1}, we know there exists a constant $\nu\in (0,1)$ depending only on $\{s,N,p,q,\lambda,\Lambda\}$, which means we shall choose $k_{1}$ and $k_{2}$ to satisfy the following conditions
    \begin{equation*}
        \frac{\gamma_{1}g\left(\frac{k_{1}}{r^{s}}\right)}{\bbint{B_{r}}g\left(\frac{u(x)}{r^{s}}\right)\differential x}\leq \nu
        \qquad \text{and} \qquad
        \frac{\gamma_{2}g\left(\frac{k_{2}}{r^{s}}\right)}{r^{s}\Tail(u_{+};B_{r})}\leq \nu,
    \end{equation*}
    where $\gamma_{1}$ and $\gamma_{2}$ are determined by Lemmas \ref{lemma2} and \ref{lemma3}. Utilizing the inequality \eqref{G,a<1}, we can choose
    \begin{equation*}
        k_{1}= r^{s}g^{-1}\left[\frac{\nu}{\gamma_{1}}\bbint{B_{r}}g\left(\frac{u(x)}{r^{s}}\right)\differential x\right]=\eta_{1}r^{s}g^{-1}\left[\bbint{B_{r}}g\left(\frac{u(x)}{r^{s}}\right)\right],
    \end{equation*}
    and
    \begin{equation*}
        k_{2}= r^{s}g^{-1}\left(\frac{\nu}{\gamma_{2}}r^{s}\Tail(u_{+};B_{r})\right)=\eta_{2}r^{s}g^{-1}\left(r^{s}\Tail(u_{+};B_{r})\right),
    \end{equation*}
    where $\eta_{1}$ and $\eta_{2}$ can be described respectively by
    \begin{equation*}
    \eta_{1}=\left(\frac{q}{p}\right)^{\frac{1}{p-1}}\left(\frac{\nu}{\gamma_{1}}\right)^{\frac{1}{q-1}}
    \qquad \text{and} \qquad
    \eta_{2}=\left(\frac{q}{p}\right)^{\frac{1}{p-1}}\left(\frac{\nu}{\gamma_{2}}\right)^{\frac{1}{q-1}}.
    \end{equation*}
    It comes to us from Lemma \ref{lemma1} that for $i=1,2$, either
    \begin{equation*}
        r^{s}g^{-1}(r^{s}\Tail(u_{-};B_{R}))+r^{s}G^{-1}G^{*}(d_{1}r^{s})>k_{i},
    \end{equation*}
    or
    \begin{equation*}
        u\geq \frac{1}{2}k_{i}.
    \end{equation*}
    Therefore, the proof is finished.
\end{proof}

\begin{thebibliography}{[a]}
\bibitem{P. Baroni2015} P. Baroni, Riesz potential estimates for a general class of quasilinear equations, Calc. Var. Partial Differ. Equations 53 (2015), 803-846.
\bibitem{bögelein2024gradientregularityspharmonicfunctions} V. B\"ogelein, F. Duzaar, N. Liao, G. Bisci and R. Servadei, Gradient regularity for $(s,p)$-harmonic functions, arXiv:2409.02012.
\bibitem{bogelein2024regularity} V. B\"ogelein, F. Duzaar, N. Liao and G. Bisci, R. Servadei, Regularity for the fractional $p$-Laplace equation, arXiv:2406.01568.
\bibitem{bonder2023global} J. F. Bonder, A. Salort, and H. Vivas, Global H\"older regularity for eigenfunctions of the fractional $g$-Laplacian, J. Math. Anal. Appl. 526.1 (2023), 127332.
\bibitem{WOS:000447961200017} L. Brasco, E. Lindgren, and A. Schikorra, Higher H\"older regularity for the fractional $p$-Laplacian in the superquadratic case, Adv. Math. 338 (2018), 782--846.
\bibitem{buryachenko2020local} K. Buryachenko and I. Skrypnik, Local continuity and Harnack’s inequality for double-phase parabolic equations, Potential Anal. 56 (2022), 137--164
\bibitem{byun2024holder} S. Byun and K. Kim, A H\"older estimate with an optimal tail for nonlocal parabolic $p$-Laplace equations, Ann. Mat. Pura Appl. (1923-) 203.1 (2024), 109--147.
\bibitem{byun2021local} S. Byun, H. Kim, and J. Ok, Local H\"older continuity for fractional nonlocal equations with general growth, Math. Ann. 387 (2022), 807--846.
\bibitem{byun2023nonlocal} S. Byun, H. Kim, and K. Song, Nonlocal Harnack inequality for fractional elliptic equations with Orlicz growth, Bull. Lond. Math. Soc. 55 (2023), 2382--2399.
\bibitem{BYUN2022110} S. Byun, J. Ok, and K. Song, H\"older regularity for weak solutions to nonlocal double phase problems, J. Math. Pures Appl. 168 (2022), 110--142,
\bibitem{2016Regularity} M. Cozzi, Regularity results and Harnack inequalities for minimizers and solutions of nonlocal problems: a unified approach via fractional De Giorgi classes, J. Funct. Anal. 272.11 (2016), 4762--4837.
\bibitem{de2019holder} C. De Filippis and G. Palatucci, H\"older regularity for nonlocal double phase equations, J. Differ. Equ. 267.1 (2019), 547--586.
\bibitem{di2016local} A. Di Castro, T. Kuusi, and G. Palatucci, Local behavior of fractional $p$-minimizers, Annales H. Poincaré C, 33.5 (2016), 1279--1299.
\bibitem{L. Diening2009} L. Diening, B. Stroffolini,and A. Verde, Everywhere regularity of functionals with $\varphi$-growth, Manuscripta Math.
129.4 (2009),449-481.
\bibitem{2022Harnack} Y. Fang and C. Zhang, Harnack inequality for the nonlocal equations with general growth, P. Roy. Soc. Edinb. A 153.5 (2023), 1479--1502.
\bibitem{garain2024higher} P. Garain and E. Lindgren, Higher H\"older regularity for the fractional $p$-Laplace equation in the subquadratic case, Math. Ann. 390 (2024), 5753--5792.
\bibitem{giusti2003direct} E. Giusti, Direct methods in the calculus of variations, Scientific World Journal, 2003.
\bibitem{iannizzotto2016global} A. Iannizzotto, S. J. Mosconi, and M. Squassina, Global H\"older regularity for the fractional $p$-Laplacian, Rev. Mat. Iberoam. 32.4 (2016), 1353--1392.
\bibitem{kassmann2011new} M. Kassmann, A new formulation of Harnack’s inequality for nonlocal operators, Comptes Rendus. Math. 349.11-12 (2011), 637--640.
\bibitem{kassmann2022nonlocal} M. Kassmann and M. Weidner, Nonlocal operators related to nonsymmetric forms \uppercase\expandafter{\romannumeral1}: H\"older estimates, arXiv:2203.07418.
\bibitem{kassmann2024nonlocal} M. Kassmann and M. Weidner, Nonlocal operators related to nonsymmetric forms \uppercase\expandafter{\romannumeral2}: Harnack inequalities, Anal. PDE 17.9 (2024), 3189-3249.
\bibitem{kassmann2024parabolicharnackinequalitynonlocal} M. Kassmann and M. Weidner, The parabolic Harnack inequality for nonlocal equations, Duke Math. J. 173.17 (2024), 3413-3451
\bibitem{kim2023supersolutionssuperharmonicfunctionsnonlocal} M. Kim and S. Lee, Supersolutions and superharmonic functions for nonlocal operators with Orlicz growth, arXiv:2311.01246.
\bibitem{liao2024nonlocalweakharnackestimates} N. Liao, Nonlocal weak Harnack estimates, arXiv:2402.11986.
\bibitem{Regularityof2017}J. Ok, Regularity of $\omega$-minimizers for a class of functionals with non-standard growth, Calc. Var. Partial
Differ. Equations 56.2 (2017), 48.
\bibitem{Glaplacemeasures2015}J. Zheng, B. Feng, Z. Zhang, Regularity of solutions to the $ G $-Laplace equation involving measures, Z. Anal. Anwend. 34.2 (2015), 165-174.
\end{thebibliography}
\end{document}